1. Basic concept of programming, data types, vectors and factors

2. Matrices and data.frames, data import/export, graphics and visualization

3. Tidy data, data wrangling

4. Today: Spatial data in R

library(raster)

ras <- raster(ncol = 100, nrow = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)

ras

## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

m <- rnorm(n = ncell(ras), mean = 0, sd = 1)

values(ras) <- m

plot(ras)

ras2 <- ras + 10

ras3 <- ras * ras2

ras4 <- ras^2

ras5 <- (ras - cellStats(ras, stat = mean)) / 
                                cellStats(ras, stat = sd)

ras.na <- ras

ras.na[ras.na > 0.5] <- NA

k <- matrix(c(  0, 0.5, 1,
              0.5, 1.5, 2,
              1.5, 6.0, 3), 
            ncol = 3, byrow = TRUE)
k

##      [,1] [,2] [,3]
## [1,]  0.0  0.5    1
## [2,]  0.5  1.5    2
## [3,]  1.5  6.0    3

ras.reclass <- reclassify(x = ras, rcl = k)

ras[10, 20]

##           
## -2.136716

ras[10, 20] <- 10

ras.stack <- stack(ras, ras2)

ras.stack

## class       : RasterStack 
## dimensions  : 100, 100, 10000, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## names       :   layer.1,   layer.2 
## min values  : -3.497505,  6.502495 
## max values  :   10.0000,   14.1605

layer1 <- subset(x = ras.stack, subset = 1)

layer1

## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer.1 
## values      : -3.497505, 10  (min, max)

data <- as.data.frame(ras.stack)

head(data)

##      layer.1   layer.2
## 1 -1.2699512  8.730049
## 2  0.5653630 10.565363
## 3 -0.3784131  9.621587
## 4 -0.6302751  9.369725
## 5 -0.1924782  9.807522
## 6  0.2749651 10.274965

elevation <- raster("Data/elevation.envi")

elevation

## class       : RasterLayer 
## dimensions  : 709, 1306, 925954  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 16.90833, 27.79167, 44.33333, 50.24167  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/corneliussenf/Dropbox/Teching/Quantitative Methods/Week-13/Data/elevation.envi 
## names       : elevation

projection(elevation)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

See for projection definitions: http://spatialreference.org/

plot(elevation)

climate <- stack("Data/climate.envi")

nlayers(climate)

## [1] 4

names(climate)

## [1] "X.Tmax."  "X.Tmin."  "X.Tmean." "X.Prec."

names(climate) <- c("Tmax", "Tmin", "Tmean", "Prec")

writeRaster(elevation, 
            filename= "Data/elev.envi",
            overwrite = TRUE)

library(sp)

coords_srs <- data.frame(x = runif(100, 17, 27), 
                         y = runif(100, 45, 50))

pts_srs <- SpatialPoints(coords = coords_srs)
pts_srs

## class       : SpatialPoints 
## features    : 100 
## extent      : 17.02223, 26.78149, 45.0296, 49.96563  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

proj <- projection(elevation)

CRS(proj)

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

pts_srs <- SpatialPoints(coords = coords_srs,
                    proj4string = CRS(proj))

plot(elevation)
plot(pts_srs, add = TRUE)

library(mapview)

mapview(climate)

df.climate <- extract(climate, pts_srs, df = TRUE)

class(df.climate)

## [1] "data.frame"

head(df.climate)

##   ID     Tmax      Tmin    Tmean     Prec
## 1  1 140.4167 32.583332 86.25000 53.66667
## 2  2 105.2500 23.500000 64.16666 64.50000
## 3  3 117.9167 30.000000 73.66666 68.58334
## 4  4 145.7500 40.000000 92.75000 52.66667
## 5  5  91.5000  9.166667 50.16667 71.58334
## 6  6 116.3333 22.916666 69.41666 64.08334

spdf.climate <- extract(climate, pts_srs, sp = TRUE)

spdf.climate

## class       : SpatialPointsDataFrame 
## features    : 100 
## extent      : 17.02223, 26.78149, 45.0296, 49.96563  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 4
## names       :             Tmax,              Tmin,            Tmean,             Prec 
## min values  : 61.9166679382324, -16.4166660308838, 22.5833339691162, 43.4166679382324 
## max values  :           166.75,  67.9166641235352,              117, 84.5833358764648

You can acess SpatialPointsDataFrame the same way as data.frame!

head(spdf.climate)

##       Tmax      Tmin    Tmean     Prec
## 1 140.4167 32.583332 86.25000 53.66667
## 2 105.2500 23.500000 64.16666 64.50000
## 3 117.9167 30.000000 73.66666 68.58334
## 4 145.7500 40.000000 92.75000 52.66667
## 5  91.5000  9.166667 50.16667 71.58334
## 6 116.3333 22.916666 69.41666 64.08334

spdf.climate$Tmax[1]

## [1] 140.4167

library(rgdal)

pts.sample <- readOGR('Data', 'pts_sample')

## OGR data source with driver: ESRI Shapefile 
## Source: "Data", layer: "pts_sample"
## with 142 features
## It has 1 fields
## Integer64 fields read as strings:  id

pts.sample

## class       : SpatialPointsDataFrame 
## features    : 142 
## extent      : 17.80209, 27.04407, 44.64419, 49.84834  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : id 
## min values  :  0 
## max values  : 99

writeOGR(spdf.climate, "Data", "pts_climate", driver = "ESRI Shapefile", overwrite_layer = TRUE)